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ABSTRACT 

In this paper we address the issue of the origin of LBV bipolar bubbles. Previous 
studies have explained the shapes of LBV nebulae, such as rj Car, by invoking the 
interaction of an isotropic fast wind with a previously deposited, slow aspherical wind 
(a "slow torus"). In this paper we focus on the opposite scenario where an aspherical 
fast wind expands into a previously deposited isotropic slow wind. Using high 
resolution hydrodynamic simulations, which include the effects of radiative cooling, we 
have completed a series of numerical experiments to test if and how aspherical fast 
winds effect wind blown bubble morphologies. Our experiments explore a variety of 
models for the latitudinal variations of fast wind flow parameters. The simulations 
demonstrate that aspherical fast winds can produce strongly bipolar outflows. In 
addition the properties of outflows recover some important aspects of LBV bubbles 
which the previous "slow torus" models can not. 

Subject headings: hydrodynamics - stars: mass loss - stars: supergiants 
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1. Introduction 



In just a few years the HST has transformed our understanding of the massive unstable stars 
know as Luminous Blue Variables (LBVs). Recent observations have revealed a number of LBVs 
or LBV candidates to be surrounded by extended aspherical outflows. The most extraordinary 



of these is the markedly bipolar nebula surrounding rj Carinae ( "the homunculus" : Hester et al. 



1991 ; Ebbets et al. 1993 ; Humphreys &: Davidson 1994 ). Other LBVs show nebulae with varying 
degrees of asphericity from elliptical (R127: |Nota et al. 1995 ) to strongly bipolar (which we define 
though the presence of an equatorial waist: HR Carinae: Nota et al. 1995| ; Weis et al. 1996 ). 



These morphologies are quite similar to what has been observed in Planetary Nebulae (PNe) 
which arise from low mass stars (Balick 1987; Manchado et al. 1996). The aspherical shapes of 
PNe have been successfully explained through a scenario termed the "Generalized Interacting 



Stellar Winds" model (GISW: [Kwok, Purton k Fitzgereld 1978| , [Kwok 1982| , |Kahn 1983| , [Frank 



et al. 1993| ; [Frank fc Mellema 199^ [Mellema fc Frank 1995[ ). In the GISW model an isotropic 
fast wind from the central star (a proto-white dwarf) expands into an aspherical (toroidal) slow 
wind ejected by the star in its previous incarnation as a Asymptotic Red Giant. High densities in 
the equatorial plane constrain the expansion of the fast wind. The expanding shock which forms 
quickly assumes an elliptical prolate geometry. If the ratio of mass density between the equator 
and pole (a parameter we call q, qp = Pel Pp) is high enough, then the elliptical bubble eventually 
develops a waist and becomes bipolar. 

The similarity of PNe and LBV nebulae has led to the suggestion that both families of objects 
are shaped in similar ways. In [Frank, Bahck, Davidson 1995[ (hereafter: FBD) a GISW model 
for T] Car was presented in which a spherical "outburst" wind expelled during the 1840 outburst 
expanded into a toroidal "pre-outburst" wind. FBD showed that the resulting bipolar outflow 



could recover both the gross morphology and kinematics of the Homunculus. Nota et al. 1995 
(hereafter NLCS) used a similar model for other LBV nebulae presenting a unified picture of the 
development of LBV outflows. More recently Garcia-segura et al. 1997[ (hereafter GLM) presented 
a model which also relied on the GISW scenario but which changed the order of importance of 
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the winds. The novel aspect of GLM's study was to include the effects of stellar rotation. Using 
the Wind Compressed Disk model of pjorkman &: Cassinelli 1992 , GLM showed that a strong 



equator to pole density contrast would likely form during the outburst when the star is close to 
the Eddington limit and rotation can deflect wind streamlines toward the equator. In their model 
it is the "post-outburst" mass loss (which was not considered in either FBD or NLCS) that acts 
as the fast wind. The post-outburst wind in GLM's model "inflates" the bipolar bubble via its 
interaction with the toroidal outburst wind. 

All these models have demonstrated the potential efficacy of the GISW scenario by recovering 
the basic shapes observed in LBV nebulae. However, by relying on the presence of a slow torus 
they are are troubling in their mutual inconsistency. Specifically the question "Where is the 
torus?" must be answered. Does the torus form during the outburst phase as in GLM or does it 
form in the pre-outburst wind as in FBD and NLCS? Without invoking binary interactions or a 
pre-existing disk left over from the stellar formation process, it may be difficult to get a strong 
toroidal density contrast in the pre-outburst environment. 

Stepping back further one can also ask if a disk is needed at all? The latter question arises 
from consideration of new HST images of r] Car orse et al. 1997| ) which reveal the disk to be 



so highly fragmented that it may be more reasonable to consider the structure to be a "skirt" 
of individual clumps of ejecta rather than a continuous feature. This point is crucial since a 
discontinuous equatorial spray of isolated bullets can not hydro dynamically constrain an isotropic 
stellar wind to form a bipolar outflow. Thus we are led to consider an alternative model to the one 
proposed by FBD, NLCS and GLM which further generalizes the GISW model by turning that 
scenario on its head. It what follows we consider the case of an aspherical fast wind interacting 
with an isotropic slow wind. We imagine a fast wind ejected with higher velocity along the poles 
than along the equator. The question we wish to answer is can such a wind, expanding into an 
isotropic environment, account for the shapes of LBV nebulae. 

There are a number of reasons for pursuing this line of investigation some of which were cited 
above as criticisms of the "classic" GISW model. More importantly, however, theoretical models 
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admit the possibility of aspherical fast winds in massive stars. Pauldrach & Puis 1990 have shown 
that a discontinuity (bistability) in mass loss and velocity occurs when the effective gravity of the 
star drops below a critical value. [Lamers fc Pauldrach 1991 used these results to demonstrate 



that stellar rotation can induce latitudinal changes in (jg// and the optical depth of the wind. 
The change in optical depth puts the polar and equatorial regions of the star on different sides 
of the bistabilty limit. A high velocity, low density wind forms at the poles, and low velocity, 
high density wind forms along the equator. It should be noted that the Wind Compressed Disk 
(WCD) model of Bjorkman & Cassinelli 1992 also produces aspherical winds since the equatorial 



focusing occurs close to the star. Thus a wind that has been shaped by the WCD mechanism, if 
it is expanding into a slower moving environment, should be considered an aspherical fast wind, 
i.e. the issue is always the velocity (and density) of previously ejected material. It is worth noting 



however that recent numerical models of the WCD mechanism ( Owocki, Cranmer fc Gayley 1996 ) 
which include non-radial line forces found inhibition of the wind compression and mass loss in the 
equator. Instead a net flow in the pole-ward direction was formed. This is, therefore, yet another 
way by which fast winds might become aspherical. Finally, and most importantly, there is direct 



evidence for asphericity in fast winds. Observations of the wind of AG Carinae (Leitherer et al 



1994) imply a pattern of densities and velocities from pole to equator much like that described in 



piiamers &: Pauldrach 199"l . Finally we note that it is worthwhile to pursue this kind of investigation 



simply because it has not been done before. The GISW model and its variations has been very 
successful in accounting for a variety of bipolar outflow phenomena ( pVIellema, Eulderink fc Icke 



1991; Blondin & Lundqvist 1993; [Frank, Balick fc Livio 1996 ). Since the effect of aspherical fast 



winds has yet to be investigated the potential of finding useful results is high which argues for a 
detailed study. 

We note that this paper represents an initial exploratory study. We are using an admittedly 
ad-hoc formalism to control the asphericity of the fast wind and we have not tuned our parameters 
to the accepted values for any particular LBV. Our purpose in this paper is to map out the broad 
consequences of including aspherical fast winds into the GISW formalism with an application to 
LBVs as a class of outflows. In future papers we will attempt to apply our results to individual 
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LBVs in an attempt to make detailed contact with observational results. 

The organization of the paper is as follows: In section II we describe the numerical method 
and initial conditions used in our simulations. In Section III we present and discuss the results of 
our simulations. In section IV we present our conclusions along with a discussion of some issues 
raised by the simulations. 



2. Computational Methods and Initial Conditions 

The gasdynamic interactions we wish to study are governed by the Euler equations with a 
'sink' term in the energy equation due to radiation losses. 

§^ + V-pu = 0, (1) 

^ + V • puu = 0, (2) 

dE p2 

+ S/-u{E+p) = -^ A(T), (3) 
ot m 

where 
and 

P = P^T (5) 
m 

In the above equations ifi is the mean mass per particle and 7 is the ratio of specific heats which 
we take to be 7 = 5/3. 

We have carried out our study using the Total Variation Diminishing (TVD) method of 
Harten (1983) as implemented by Ryu et al (1995) . The TVD code is an explicit method for 



solving the Euler equations which achieves second order accuracy by finding approximate solutions 
to the Riemann problem at grid boundaries while remaining non-oscillatory at shocks through 
the application of a lower order monotone scheme. For this study we used the code configured in 
cylindrical coordinates (r, z). 
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The cooling term A(T) was calculated via look-up tables for A(T) taken from the coronal 



cooling curve of Dalgarno & McCray 1972. In the TVD code the cooling is applied in between 



hydro timesteps via an integration of the thermal energy (£"4) equation, 

If the cooling time scale becomes smaller than the dynamical time scale, the cooling term 
becomes a "stiff' source term. There are several ways to the handle stiff source terms (see, e.g., 
LeVeque 1997| ). First, the cooling time scale as well as the Courant condition is considered to 



determine the timestep. Then, in some situations, the timestep is governed by the cooing time 
scale and becomes uncomfortably too small. Second, Strang's operator splitting is employed, 
where the cooling is computed by multiple steps with it own timestep during one hydrodynamic 
timestep. Third, the solution takes the form 

Et^+' = Et^e^p(^-^At^ (7) 

where the superscript n refers to the time index. In this method, the cooling is computed in 
a single step. It produces good results even though the cooling time scale is smaller than the 
dynamical time scale. When cooling is small, the solution approaches the conventional form 

£^"+1 = Et"" - E'^M. (8) 

The third method has been used in our simulations. The term in the exponential can be 
expressed as At/At^ where At^ is the local cooling timestep of the gas At^ = E^'/E^'. Although 
the method works well for small cooling timestep, for safety we have used 

At = min[At/i,1.5Atc], (9) 

where Ath is the dynamical timestep set by the Courant Condition. 

Care must be taken in the application of cooling near contact discontinuities (CDs) where 
in typical wind blown bubble temperatures and densities can vary by more than an order of 
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magnitude going from low to high values of density and high to low values of temperature. In 
simulations with less than perfect resolution the CD will be smeared out across a number of cells. 
In these intermediate zones which have both high density and high temperature the cooling term 
can be anomalously large. This leads to the removal of prodigious amounts of energy from the 
system at just a few zones. To avoid these problems we have tracked the CD by following the 
advection of a passive tracer i; via Eq. @, i.e. 

^ + V • Vu = 0. (10) 

This equation has been solved with the TVD code. Using this tracer we can distinguish between 
shocks and CDs. If the fast stellar wind has ip = 1 and the slow material initially filling the grid 
has ip = then a CD can always be marked as a location with 0.1 < "0 < 0.9. When a CD is 
encountered we replace the anomalous cooling value with one obtained from an average of the 8 
nearest neighbor cells which do not contain the CD. Tests of the code show that this formalism 
allows us to reproduce analytically derived shock speeds ( Koo Sz McKee 1992| ) for spherical 



radiative bubbles in a variety of power-law, p cc r^, environments. 

Our numerical experiments are designed to study the evolution of wind-blown bubbles driven 
by aspherical fast winds. The environment is always assumed to be characteristic of a previously 
deposited spherically symmetric "pre-outburst" wind which we denote as wind with mass loss 
rate Mq and velocity Vq. Thus po = Mq/ Airr'^VQ. For the driving "fast" or "outburst" wind, which 
we denote as wind 1, we need a formalism for setting the latitudinal (6) variations in the wind 
properties, i.e. Mi = Mi{6) and Vi = Vi{9). We note that since we wish to drive prolate bipolar 
bubbles we always assume that the velocity at the poles is larger that at the equator 

We have chosen to explore 3 models for the pole to equator variation in wind parameters. 
Each model is based on the assumption of a different quantity remaining constant across the face 
of the star. Our three models are: 



1. Constant Momentum Input 
• n = MiFl = Const 
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2. Constant Energy Input 

• E = \MiVi^ = Const 

3. Constant density 

• pi = Const 

If we choose our fiducial values for the density and velocity at the equator (pie, Vie), then the 
variation of these quantities can be expressed as pi{9) = pieg{0), and Vi{9) = Vief{6). We can 
then write g{9) oc f{9)^. For case 1 where 11 = Const = He, we have 

n = Mi{e)Vi{e) = (4^i?,Vi(^)Vi(e))^i(0) = (47ri?,Vie^ie)^ie. (n) 

Thus 



9{G) = TTZ^- (12) 



1 

m 

For case 1 n = — 2. Similarly, for cases 2 and 3 one finds n = — 3 and respectively. 



For f[9) we have chosen an ad hoc function which produces a smooth variation in Vi and p\ 
from equator to pole. Written in terms of velocity we have 

exp(-2Scos(^)2) - 



Vi{e) = Vi 



le 



l-A- 



(13) 



exp {-2B) - 1 

where constant B determines the shape of the equator to pole variability. The constant A is related 
to the magnitude of the contrast A = \ — q^. Note in what follows we shall write qx = X(,/Xp and 
X is either the velocity or the density. Note that Vp = 14/(1 — A)soVp = Ve/qv We take A <1, 
so that together with Eq. (p^) we have q^ = Ve/Vp < 1 and qp = Pe/ Pp > 1- These relations express 
our assumptions about the velocity variation across the face of the star and its consequences for 
the density distribution given our different models. In Fig. 1 we show a plot of the fast wind V{0) 
and p{6) for each of our three cases. 

Using the formalism described above we have carried out three sets of numerical experiments 
varying the assumptions about the mass loss rates in the successive winds between each set. 
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Within the first two sets we performed a triplet of simulations with the 9 dependence of Mi and 
Vi corresponding to the three cases discussed above. In set 1 we examined the interaction between 
two winds of with equal mass loss rates. These simulations are performed to give us a baseline 
on the gas-dynamic flow pattern. In the second set of experiments we examined the interaction 
between two winds where the previously ejected "slow" pre-outburst wind (Mq) had a lower mass 
loss rate than the outburst wind. We examine this case because stellar evolution models show 
the mass loss rate increases during an LBV eruption ( Langer et al. 1994| ). In the third set of 



experiments we examined the interaction of three winds where Mq = M2 < Mi and Vq < Vi < V2. 
In these simulations only the outburst wind is considered to be aspherical. We examine this case 
to explore the effect of the post outburst wind on the nebular morphology. We note that the 
velocity difference between pre and post outburst winds may not be physically realized. As we 
shall discuss in a section 3.3 this assumption should not effect our results significantly. The initial 
conditions for each of our 9 simulations including the index n of the density function {i.e. Eq. [l^) 
are given in table 1. 

The mass loss rates use in our simulations range from M = 10~^ M0 yr~"'^to M = 10~^ 
M© yr""*^. Our wind velocities vary between 100 km s^^and 1400 km s^^. These values are 



representative of what is observed in LBVs (Leitherer et al. 1994) with the highest velocities seen 



in rj Car ( Ebbets et al. 1997| ). We note that in giant outbursts the mass loss rates may increase 



by an order of magnitude or more over what is used in this study ( [Humphreys &i Davidson 1994| ). 
We have chosen not to use larger values of the mass loss rate based on the short cooling timescales 
which occur in high density winds. Short timescales would force us to use lower resolution 
simulations given constrains on computational time. We decided that achieving higher resolution 
(256 X 256) was a more important goal since, based on our results, we will be able to anticipate 
the enhanced cooling effects of higher mass loss rates. We note that the simulations were all run 
until the outer shock ran off the grid. Thus the timescale for each simulation is essentially the 
dynamical age of the bubble for i? ~ 5 x lO^'' cm. 
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3. Results 

3.1. 2 Wind Models with Equal M 

Fig. 2 shows /3, T, V ^ and P for run A after t = 800 yrs of evolution. Run A has an outburst 
stellar wind whose momentum flux n is constant across the face of the star. A number of features 
of the simulation are noteworthy but it is first worthwhile to identify the main features of the wind 
blown bubble. The density map shows a bright rim which corresponds to the shell of swept-up 
pre-outburst wind material. It is important to distinguish the composition of the material in the 
shell (i.e. which wind it originates from) since its high density will make it most luminous and, 
therefore, the defining feature in a real nebula. Fig. 2 shows that the shell develops significant 
asphericity due to the aspherical driving force of the outburst wind. For later comparison we 
introduce a quantity we call the ellipticity defined as the ratio of the distance from the star to 
the outer shock in the equator and the pole, e = Re/Rp- e ~ 0.68 in this simulation. Notice also 
that the bubble has developed a "waist" which is the observational signature of a bipolar rather 
than an elliptical configuration. The density map for Run A is also shown in Fig. 4 with the 
appropriate reflections to create a slice of the full bubble. 

Detailed inspection of the p and T maps show that the shell is thin and cold (T ~ W^K) at 
all latitudes. The narrowness of the shell is due to the enhanced compression which accompanies 
radiative losses behind the outer shock. The role of cooling is important and is determined by 
the relation between the age of the bubble and the cooling timescale. Given a shock speed Vs 
and a preshock density of ppre the cooling timescale can be defined as tc = CVs'^/ppre where 
C = 6 X 10^^^ g cm~^ s"^ (Kahn 1976). The outer shock speed in this simulation is on the order of 
200 km/s. Taking pp^e from the mass conservation in the pre-outburst wind gives tc = 150y which 
is short compared to the age of the bubble in Fig. 2. 

The situation is quite different behind the "reverse "or "inner" shock. At low latitudes 
the speed of the outburst wind is relatively small (150 km/s) and the density is relatively high 
producing a short cooling time. The reverse situation occurs at the poles where the high wind 
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velocity (750 km/s) and low pre-shock density yield cooling times which are of the same order or 
longer than the age of the bubble. This produces a cap of high temperature gas at the poles. High 
thermal pressure enhances the preferential expansion at the poles which leads to a bipolar bubble. 
We note however that even when the cooling at the poles is strong enough to drain away all shock 
thermalized energy the bubble which is produced is still bipolar. 

The density maps for run B and C are presented in Fig. 4. Run B corresponds to the case 
where the energy input in the wind, Ei, is constant across the face of the star. The morphology 
of the bubble in run B is similar to that in run A. The principle difference being the aspect ratio 
of the bubble. This can be attributed to the differing latitudinal density profiles for the two runs. 



From Eq. (12) the ITi = Const simulation (run A) has p oc f{9)^^. The Ei = Const simulation 
(run B) yields p oc f{9)~^. For run C the density in the wind is constant and is set by the value 
in the equator where Vi{9) has its minimum. Since p oc M /V via the continuity equation the 
outburst wind density in run C is almost as high as that in the pre-outburst wind. This produces 
a strong ram pressure flux ratio between the two flows. The initial shape of the outburst wind is 
strongly imprinted on the developing bubble. Thus run C produces most bipolar shell of all three 
simulations examined in this first set of experiments. 

The results presented above show that bipolar bubbles can form when driven by an aspherical 
fast wind. The differences between the three models also demonstrates that the degree of bipolarity 
does depend on relative densities between the aspherical fast (outburst) wind and the spherically 
symmetric slow pre-outburst wind. 



3.2. 2 Wind Models with Mi > Mq 



As [Langer et al. 1994 have demonstrated with stellar evolution codes, the LBV outburst phase 



is likely to involve an increase in mass loss over the pre-outburst wind. Thus we have the case of a 
"heavy" wind expanding into a light one. In this situation the pre-existing circumstellar material 
can not constrain the outflow. The wind will not be significantly decelerated by the environment 
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until the swept-up mass is comparable to the mass ejected by the wind. Given that the outburst 
mass loss rate can be an order of magnitude or more greater than the pre-outburst value, the 
outburst wind can expand almost ballistically for some time. If, in addition, the outburst wind 
is aspherical, then the bubble which is created will retain much of the asphericity imprinted on 
the wind at the stellar surface. In this scenario we expect that one gets out something similar to 
what is put in via the initial density and velocity distributions for the outburst wind. To test this 
expectation we have run a second set of experiments consisting of three simulations. Runs D, E 
and F have Hi, Ei and pi held constant respectively just as in the simulations described in the 
last section. In these simulations however Mi = IOOMq. Note that for the sake of consistency with 
the previous experiments we have used the same velocity in the pre-outburst wind as was used in 



runs A through C. The results of Langer et al. 1994 show that the velocity in the pre-outburst 



wind should be higher than that in the outburst wind. We do not believe such a difference will 
produce a profound change in our results. A high velocity pre-outburst wind will interact with 
the ISM producing shocked gas which backflows to the boundary of the pre-outburst wind. Such 
a high pressure region will still yield an isotropic force opposing the expansion of the outburst 
wind. Thus while the details of the flow may change, the global properties most likely will not be 
affected. 

Fig. 3. shows a snapshot of the evolution of run D after t = 240yrs. As expected the bubble 
formed appears strongly bipolar with an ellipticity of e = 0.4. While such a result was expected 
from intuitive arguments a ballistic expansion model does not tell the whole story. Note first that 
run D has reached the same size {Re = 5.1 x 10^^ cm) as run A in about a third of the time giving 
an expansion velocity along the pole of Vg ~ 660 km/s. This is less than the outburst wind velocity 
at the poles Vi = 750 km/s. Ambient material which passes through shock can be expected to 
reach high temperatures of the order of T ~ 10^ K where cooling is relatively ineffective. The 
temperature map clearly shows high T gas at the top of the bubble. Consideration of the passive 
tracer shows that this is ambient material which has been accelerated. The temperature map also 
shows lower temperature gas at smaller radii. This is outburst wind material which has passed 
through the inner shock which can be seen at Re = 4.2 x 10^'^ cm. Consideration of the position 
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of the shock compared with the position expected for baUistic winds yields a inner shock velocity 
of Vg ~ 200 km s~^, consistent with the temperature in this region (measured to T ^ 10^ K). The 
presence of this shock indicates that the outburst wind gas parcels arc being decelerated by the 
swept-up ambient material. This deceleration occurs because of the low outburst wind density 
along the poles as can see from the density map. Recall that we are using equator to pole velocity 
contrast of qy = 1/5. This produces a equator to pole density contrast of qp = l/qy^ = 25. 

Runs E and F show similar morphologies though the p = Const case again produces the most 
bipolar configuration. We note that for Run E and F e = 0.42 and 0.36 respectively. In Fig. 4 
we present logio P maps for runs A through F. In each map we have reflected the computational 
space about the symmetry axis (r = 0) and symmetry plane {z = 0) to show the full bubble cross 
section and facilitate comparison. In spite of the differences between the runs, Fig. 4 demonstrates 
the principle conclusion of our first two sets of experiments: an aspherical stellar wind can drive 
an aspherical bubble. 



3.3. 3 Wind Models 

In this section we present the results of our final set of numerical experiments where we take 
a step closer to reality by invoking three wind models. In runs G, H, and I we have performed 
simulations with conditions in the pre-outburst and outburst winds similar to run D. In these 
runs however the outburst wind lasts only 30 years after which a "post-outburst" wind is driven 
into the grid. The characteristics of the post-outburst wind were meant to mimic the conditions 
currently observed in r] Carinae in the sense that we used a relatively low mass rate and high 
velocity for the wind, i.e. M2 < Mi and V2 > Vi. Recall that in the simulations of GLM it was 
the fast post-outburst wind which was responsible for driving the bipolar rj Car bubble. In these 
simulations we are interested in what effect this post-outburst wind will have on a bipolar bubble 
created by a previously ejected dense aspherical wind. We note that the post-outburst wind in all 
of these simulations is spherical and will produce an isotropic driving force. 
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We have run three simulations each of which utihzed the n = Const formahsm. In these 
simulations it was q^, the equator to pole velocity contrast in the outburst wind, which was varied. 
The purpose of this strategy was to understand how the additional driving force provided by the 
post-outburst wind would effect bubbles with different degrees of bipolarity. 

Fig. 5 shows a snapshot of the evolution of run H after t = 250yrs. This simulation has 
an outburst wind with = 1/7 and Qp = 49. Unlike the last set of experiments here the gas 
ejected during the outburst in these models occupies a thin shell. The reasons for this are twofold. 
First, the outburst had a short intrinsic lifetime (which we denote as ti) so its intrinsic width 
is 6R = Viti- Second, the outburst material has been decelerated and compressed from the 
outside by the pre-outburst wind and compressed and accelerated from the inside by the action 
of succeeding post-outburst wind. Note that the post-outburst wind has itself been decelerated 
through the action of a strong inner shock. The temperature behind this shock at the poles is 
quite high (T ~ 2. x 10^ K) reflecting the high initial velocity of the post-outburst wind. The inner 
shock has a highly aspherical configuration. The radially streaming post-outburst wind encounters 
this shock at oblique angles and is refracted towards the symmetry axis ( [Frank fc Mellema 1996] ). 
This produces strong shearing flows in the region behind the inner shock which appear to give rise 
to instabilities in the interface separating the shocked post-outburst and outburst flows. 

In Fig. 6 we present log^^g P maps for all three simulations from this set of experiments. 
Again it is clear that strong bipolar morphologies develop without the need for a slow-moving 
disk. Comparison between the simulations shows that decreasing produces stronger bipolar 
morphologies. It is worth noting that if the expansion of the bubble were ballistic then we would 
expect e q^ since Re = Viet and Rp = Vipt. This is not the case. Thus these bubbles show the 
result of significant hydrodynamic shaping. Some part of the change in the shape results from the 
deceleration of the outburst wind via the previously ejected material. The post-outburst wind 
however also contributes by accelerating the outburst wind material. As the bubble evolves this 
acceleration will have its greatest effect near the equator where the outburst wind has been most 
strongly decelerated. Thus the action of the post-outburst wind will be to drive the bubble towards 
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a more spherical configuration as system evolves. We note however that the post-outburst might 
also be aspherical. As was noted above [Leitherer et al. 1994 have observed two components in the 
wind of AG Car. An strong aspherical post outburst wind would then increase the elongation of 
the bipolar lobes. Further observational study of the characteristics of LBV winds are warranted 
to determine the latitudinal density and velocity variations. 



4. Discussion and Conclusions 

The results of our simulations demonstrate that bipolar wind blown bubbles can result purely 
from the action of an aspherical fast wind. In previous studies of LBV bubbles (FBD, NLCS, 
GLM) it has been assumed that a slow moving torus or disk of gas was a necessary precondition 
for the development a bipolar bubble. Our results indicate that the properties of LBV bubbles 
may not require such a torus to form either before (FBD, NLCS) or during (GLM) the outburst. 

Our scenario has a number of attractive features. First it is both observationally and 
theoretically motivated. From the observational side there is already evidence that LBV winds 
can take on aspherical velocity and density distributions. From the theoretical side the theory of 



Lamers & Pauldrach 1991 have demonstrated that "bistable" winds are possible around massive 
hot stars. In addition, the diversity of shapes of LBV bubbles (NLCS) may be difficult to achieve 
with pre-exiting disk models. The models presented here can recover the diversity of shapes 
simply by changing though it is certainly true that this also begs the question of what drives 
the velocity contrast in the fast wind. 

Currently it is unclear what form the mass distribution takes in the lobes of LBV bubbles. 
One critical test of the different models for LBV nebula shaping would entail comparison of 
latitudinal variations of mass. If the caps of a bipolar lobe have densities that are comparable to 
that in the lobe's flanks it would present difficulties for the slow torus models. In a bipolar bubble 
resulting from a spherical fast wind driving into a slow torus the caps of the bubble should be the 
location of the lowest density. For aspherical fast winds however high or equal density in the poles 
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poses no significant problem since we are then seeing a signature of the latitudinal dependence of 
the fast wind density (consider Runs C and F in Fig. 3). We wish to note also that [Currie et al. 



1996b and Morse et al. 1997| have found that the shape of 77 Car is best matched by by a geometry 



which can be described as a "double flask" rather than an two oscullating spheres. Based on 
comparison of published results the models presented here seem to do a better job of recovering 
such a shape than either FBD, NLCS or GLM. 

It is noteworthy however that the scenario presented here would not produce the equatorial 
"skirt" seen surrounding the waist of the homunculus in r] Car. The presence of that feature 
is what motivated the original application of the GISW slow torus models. Within the current 
formulation of the aspherical fast wind model there is not a likely means of producing such a 
feature. A few points are worth noting however. The equatorial skirt is not a continuous or even 
axisymmetric feature. Thus whatever its origin it is hard to imagine that it can be the agent which 
constricts a spherical fast wind and produces the bipolar bubble. In addition there are a number of 
"spikes" extending beyond, but connecting with the homunculus that have been reported to have 
velocities in excess of 1000 km/s ( [Meaburn et al. 1996| ). They yield dynamical timescales < to 
where to is the time since the outburst. Thus it is possible that the equatorial skirt is actually a 
spray of material which ejected at some point after the outburst of 1849 and which was decelerated 
by its passage through the dense shell of the outburst wind. The non-axisymmetric distribution of 
the skirt may then be a consequence of the intrinsic pattern of ejecta or of impulsive instabilities 
which will occur when the ejecta is driven through the outburst wind. 

Regardless of the answer to this issue the results presented here show that there are two 
very different scenarios for the formation of LBV bubbles. Either they form via the interaction 
of a spherical fast wind driving into an aspherical slow wind (a slow torus) or they form from 
an aspherical fast wind driving into an isotropic pre-existing environment. This embarrassment 
of riches can be eventually be dealt with by comparing the latitudinal distributions of mass and 
momentum observed in real bipolar LBVs with what is predicted for the various models. Such 
an approach was used successfully by [Masson &: Chernin 1992 in evaluating different models of 
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molecular outflow formation in YSOs. This project is currently in progress. 
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Table 1: Initial Conditions For Runs A - H 



run 


Mo 




Vb 


Ml 






M2 


V2 


Qv 


C( 


A 


1 X 10" 


-4 


100 


1 X 10- 


-4 


150 


NA 


NA 


0.2 


-2 


B 


1 X 10" 


-4 


100 


1 X 10" 


-4 


150 


NA 


NA 


0.2 


-3 


C 


1 X 10" 


-4 


100 


1 X 10- 


-4 


150 


NA 


NA 


0.2 





D 


1 X 10- 


-6 


100 


1 X 10- 


-4 


150 


NA 


NA 


0.2 


-2 


E 


1 X 10- 


-6 


100 


1 X 10- 


-4 


150 


NA 


NA 


0.2 


-3 


F 


1 X 10- 


-6 


100 


1 X 10- 


-4 


150 


NA 


NA 


0.2 





G 


1 X 10" 


-6 


100 


1 X 10- 


-4 


150 


1 X 10"^ 


1400 


0.3 


-2 


H 


1 X 10- 


'6 


100 


1 X 10- 


-4 


150 


1 X 10-6 


1400 


0.14 


-2 


I 


1 X 10- 


-6 


100 


1 X 10- 


-4 


150 


1 X 10-6 


1400 


0.1 


-2 



Con. Index 
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FIGURE CAPTIONS 

Fig. 1 Initial conditions for aspherical fast stellar wind. The variation of stellar wind velocity 
and density with latitude are shown for a model with A = 0.1 and B = 2. Note the equator 
is at = 90°. The velocity and density at the equator are set to 1. The density plot shows 
the three different cases considered in text. Solid line: momentum input H = const across 
face of star. Dashed line: energy input E = const. Dashed-dot line: density p = const. 

Fig. 2 Run A after SOOyrs. Shown are greyscale maps of log^o number density, 

temperature, velocity (with vectors superposed), and pressure. Note that dark scales 
correspond to high values of the variables. The range of variables are as follows: 

.5 < logio(n/cm-3) < 3.7, 3.9 < \og^Q{T/K) < 6.8, 6.1 < log^^iV /cms~^) < 7.8, and 
-11.3 < log^Q{P/ dynes cm~'^) < -7.8 

Fig. 3 Run D after 24Qyrs. The maps and parameter ranges are the same as in Fig. 2 

Fig. 4 logio grayscale density maps of runs A through F. Top Left: Run A. Top Center: Run B. 
Top Right: Run C. Bottom Left: Run D. Bottom Center: Run E. Bottom Right: Run F. 
Note that here that dark scales correspond to low values of the variables. Runs A and B are 
shown at i ~ 800?/rs. Runs C through F are shown sX t ^ 250yrs. Each image is 1.2x10^® 
cm square. 

Fig. 5 Run G after 250yrs. Shown are greyscale maps of logj^g number density, 

temperature, velocity (with vectors superposed), and pressure. Note that dark scales 
correspond to high values of the variables. The range of variables are as follows: 
.5 < logio(n/cm-3) < 3.1, 3.9 < logio{T/K) < 7.5, 6.5 < logio(Wcm s'^) < 8.1, and 
-11.3 < logiQ{P/dynes cm~^) < -7.7 

Fig. 6 logio grayscale density maps of runs G through L Left: Run G. Center: Run H. Right: 
Run L Note that here that dark scales correspond to low values of the variables. Each image 
is 1.2x10^^ cm square. 
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